# Replication archive for 
# Coppock, Alexander, Donald P. Green, and Ethan Porter. 2022. 
# Does Digital Advertising Affect Vote Choice? Evidence from a Randomized Field Experiment 
# Research & Politics. 

rm(list = ls())
# setwd() to replication archive

library(DeclareDesign)
library(tidyverse)

precinct_level <- read_rds("CGP_2022_precinct_level.rds")
zip_code_level <- read_rds("CGP_2022_zip_code_level.rds")


fl_blocked_zips <-
  zip_code_level %>% 
  transmute(vb_tsmart_zip, block_id)

fl_random_assignment_function <-
  function(data) {
    new_randomization <-
      fl_blocked_zips %>%
      mutate(
        Zsim = block_ra(block_id),
        treatmentsim = as.character(block_ra(
          Zsim,
          conditions = c("video_1", "video_2")
        )),
        treatmentsim = ifelse(Zsim == 1, treatmentsim, "control")
      )
    left_join(data, select(new_randomization, vb_tsmart_zip, Zsim, treatmentsim), by = "vb_tsmart_zip")
  }

design <-
  declare_model(data = precinct_level) +
  declare_assignment(handler = fl_random_assignment_function) +
  # sharp null hypothesis of no effect
  declare_estimator(dem_two_party_share_2018 ~ Zsim, clusters = vb_tsmart_zip, label = "DIMbinary") +
  declare_estimator(
    dem_two_party_share_2018 ~ Zsim +
      dem_two_party_share_2016_nona + vote_2016_missing +
      dem_two_party_share_2014_nona + vote_2014_missing +
      dem_two_party_share_2012_nona + vote_2012_missing,
    clusters = vb_tsmart_zip,
    label = "OLSbinary"
  ) +
  declare_estimator(dem_two_party_share_2018 ~ treatmentsim,
                    clusters = vb_tsmart_zip,
                    term = c("treatmentsimvideo_1", "treatmentsimvideo_2"),
                    label = "DIMthreearm") +
  declare_estimator(
    dem_two_party_share_2018 ~ treatmentsim +
      dem_two_party_share_2016_nona + vote_2016_missing +
      dem_two_party_share_2014_nona + vote_2014_missing +
      dem_two_party_share_2012_nona + vote_2012_missing,
    clusters = vb_tsmart_zip,
    term = c("treatmentsimvideo_1", "treatmentsimvideo_2"),
    label = "OLSthreearm"
  )

# simulations <- simulate_design(design, sims = 2000)
simulations <- simulate_design(design, sims = 500)

# Original to check

# term      estimate   std.error   statistic   p.value    conf.low  conf.high       df
# 1                Z  0.0021144589 0.022487296  0.09402904 0.9252763 -0.04250567 0.04673459 98.91791
# 2                Z -0.0004433233 0.008475969 -0.05230355 0.9584026 -0.01728204 0.01639539 90.10811
# 3 treatmentvideo_1 -0.0205415171 0.024977940 -0.82238637 0.4147184 -0.07069731 0.02961427 50.56461
# 4 treatmentvideo_2  0.0255149975 0.029118158  0.87625725 0.3851601 -0.03299815 0.08402814 49.06573
# 5 treatmentvideo_1 -0.0016193103 0.010860872 -0.14909580 0.8820669 -0.02342421 0.02018559 50.92376
# 6 treatmentvideo_2  0.0006758476 0.009605703  0.07035899 0.9441904 -0.01862093 0.01997262 49.67228
# outcome
# 1 dem_two_party_share_2018
# 2 dem_two_party_share_2018
# 3 dem_two_party_share_2018
# 4 dem_two_party_share_2018
# 5 dem_two_party_share_2018
# 6 dem_two_party_share_2018

observed_estimates <- 
  get_estimates(design, data = rename(precinct_level, Zsim = Z, treatmentsim = treatment))


observed_estimates <- 
  observed_estimates %>% 
  transmute(estimator, term, estimate_obs = estimate)


p_values_df <-
  simulations %>% 
  left_join(observed_estimates) %>% 
  group_by(estimator, term) %>% 
  summarise(
    estimate_obs = unique(estimate_obs),
    p.upper = mean(estimate >= estimate_obs),
    p.twotail = mean(abs(estimate) >= abs(estimate_obs))
  )

p_values_df
